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Abstract 

The Direct Simulation Monte Carlo (DSMC) method is applied in this paper to the 
study of rarefied, hypersonic, reentry flows. The assumptions and simplifications involved 
with the treatment of ionization, free electrons and the electric field are investigated. A new 
method is presented for the calculation of the electric field and handling of charged particles 
with DSMC. In addition, a two-step model for electron impact ionization is implemented. 
The flowfield representing a 10 km/sec shock at an altitude of 65 km is calculated. The 
effects of the new modeling techniques on the calculation results are presented and discussed. 
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electric charge 
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electric field 
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activation energy 
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E c 

energy of colliding particles 

Eoo 

ionization energy 

k 

Boltzmann constant 

k(T) 

chemical rate coefficient 

m 

mass 

h 

number density 

n 

principal quantum number 

AT 

nitrogen atom 

O 

oxygen atom 

Pi 

energy of excited state 

Pr 

steric factor 

<1 

electric charge 

R 

relaxation collision number 

Ry 

Rydberg constant 

t 

time 

T 

temperature 

V 

particle velocity 

X 

distance coordinate 

c 

internal degrees of freedom 

V 

collision frequency 

CT 

reaction cross-section 

Subscripts 

e 

value for electrons 

h 

value for heavy particles 

i 

value for ions 

s 

value for particle s 
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Superscript 


* particle in excited state 

Introduction 

The Direct Simulation Monte Carlo (DSMC) method of Bird[l] has been used for more 
than two decades in the numerical simulation of physical fluid flows. The method has 
been developed to the extent that it can be successfully applied to typical engineering flows 
involving shocks in air with thermal and chemical nonequilibrium. Recent work has been 
focused on extending the method to treat flows with ionization and radiation[2,3]. This level 
of complexity is necessary to accurately predict the reentry flow fields for such vehicles as the 
Aeroassist Flight Experiment (AFE) and Aeroassisted Orbital Transfer Vehicles (AOTV). 

Certain difficulties are encountered when using DSMC with flows involving ionization. 
The presence of ions and electrons in the flow leads to electric field effects such as 
charged particle acceleration and ambipolar diffusion. Free electrons in the flow have 
extremely high velocities and collision rates compared with the heavy particles. This implies 
considerably more computational time to accurately simulate these parameters. However, 
the determination of the free electron properties is particularly important because the 
radiation environment and level of electron impact ionization are dependent on the energy 
of the free electrons. 

The method which has been employed by Bird [2] in previous simulations of flow with 
ionization involves some approximations which may introduce errors in the determination of 
free electron properties. In cases which require an accurate determination of these properties, 
an alternate method is proposed which involves the calculation of the ambipolar diffusion 
induced electric field and its influence on the motion of the charged particles. An alternative 
approach to the calculation of electron impact ionization is also introduced. This method 
involves the assumption that the reaction proceeds in two steps: electron impact excitation 
followed by ionization from the excited state. 
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The DSMC Method 


Statistical methods are used in the DSMC formulation to obtain a simulation of the 
flow of low density gases. The simultaneous computation of the trajectories of thousands of 
simulated molecules are performed in physical space. Collisions between the particles and 
interactions of the particles with the boundaries of the computational region are modeled. 
The flow properties are determined from the averages of the particle properties after a very 
large number of such interactions have taken place. A cell structure in the computational 
space is used to determine potential collision partners and for sampling of flow properties. 
The time parameter in the flow may be identified with real time and is advanced according 
to the collision frequency appropriate for the flow. An important assumption in the DSMC 
method is that the molecular motion and intermolecular collisions may be uncoupled over 
the small time step, At, used to advance the calculation. Thus, collisions between particles 
proceed until the time is advanced in the flow by At, then the particles are moved the 
distance determined by their current velocities and At. This has been shownfl] to be a 
valid assumption when the time step is less than the local average collision time 

At < 1/v (1) 

where v is the local average collision frequency. The position coordinates, velocity compo- 
nents and internal state of each molecule are stored in the computer and are modified with 
time as the molecules are followed through representative collisions and boundary interac- 
tions. In addition to changes in translational and internal energy, collisions may result in 
chemical reactions or ionization. Particles enter or exit the flow at computational boundaries 
representing free stream, a vacuum, or a known flow solution. Advantage may be taken of 
flow symmetries to reduce the number of dimensions of the grid and the number of position 
coordinates that need to be stored, but the collisions are calculated as three-dimensional 
phenomena. Although the flow is always unsteady, boundary conditions may be such that 
a steady flow is obtained as the large time state of the unsteady flow. 
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Models of the physical processes in the flow are required. The Variable Hard Sphere 
(VHS) model is used for the molecular interaction potential. The Borgnakke and Larsen 
phenomenological model is employed in the rotational and vibrational internal energy calcu- 
lations. The implementation of these models is detailed in reference 1. The nonequilibrium 
chemical reactions proceed on the basis of collision energy dependent steric factors or reac- 
tion cross-sections. Because data are not generally available for the reaction cross-sections, 
a form of the collision theory that is consistent with the VHS model[4] is used to convert 
temperature dependent rate constant data to collision energy dependent steric factors. The 
reaction rates are given in reference 5. 

For the calculations presented in this paper, a one dimensional DSMC stagnation 
streamline program is used. This program can be used to calculate a standing shock wave 
flowfield by adjusting the boundary conditions. Initially, a stagnation streamline calculation 
is performed where the boundary at one end is the undisturbed freestream while the other 
boundary is a stationary wall. When the shock reaches the desired position, molecular 
removal downstream of the shock commences. The molecular removal is conducted in such 
a way that mass, momentum and energy are conserved. It has been shown that this is 
achieved if the molecules are removed with a probability proportional to the square of 
their velocity component normal to the stream[6]. For the standing shock wave simulation, 
molecule removal is immediately adjacent to the wall and the remainder of the flow is exactly 
one dimensional. 


Modeling of Plasmas 

Several difficulties arise when the modeling of plasmas is attempted with DSMC. 
In general, the DSMC programs model only binary collisions in the flow. This is valid 
when three body collisions and multi-body charged particle interactions are insignificant in 
comparison with the two-body collisions in the flow. For the hypersonic reentry flows of 
interest in this paper, multi-body charged particle collisions may be ignored if the ionization 
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is less than 3-4 %. Charged particles are thus a minor species in the flow. Because each 
simulated particle represents an extremely large number of real particles, flow fluctuations 
in the simulation are many orders of magnitude larger than those in the real gas. These 
fluctuations in a charged gas could result in the prediction of an electric field which is far 
stronger than that in the real field. Also, the electron velocity and collision frequency are 
much higher than that for heavy particles. This implies a much smaller computational time 
step. 

In early simulations [2], a simplified procedure for the modeling of plasmas was adopted. 
To prevent the diffusion of electrons out of the flow, each electron is associated with an ion 
at the time of its formation. Collisions are calculated for the electrons as for the heavy 
particles. However, they are moved in relation to the position of their associated ion 
rather than according to their velocities. This assures charge neutrality in the flow and 
eliminates the need to calculate complex electron trajectories. The time step used in the 
simulation is the one associated with molecular collisions rather than with the electron 
collision frequency. Because the movement of the electrons relative to the ions is restricted, 
the explicit evaluation of the electric field is not required. Experimental studies indicate that 
at high altitudes (above 70 km) the effects of ambipolar diffusion are significant and should 
not be ignored[7]. Therefore, an iterative solution to account for this effect is used[6]. After 
obtaining a flow field solution, the electric field is estimated from a form of the Langmuir 
and Tonks[8] equation. 

E = (kT e / e)d(ln(n e ))l dx (2) 

This electric field is then considered in calculating the movement of the ions and a new 
flowfield solution is obtained. The calculation procedure is repeated until the solution is 
converged. 

A major problem with the above method is that the position of the electrons is 
artificially constrained (each electron associated with a specific ion). In the real flow, 
electrons from the higher temperature regions will diffuse into the lower temperature regions 
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and vice versa. The use of the Langmuir and Tonks equation is questionable because it is 
a continuum formulation. The dependence of the equation on the gradient of the electron 
number density magnifies the errors associated with the statistical sampling of a minor 
species. Equation (2), which is the form of the equation used in the previous DSMC 
simulations, can be derived from the (continuum) momentum equations for the ions and 
electrons. Its derivation requires the following assumptions; slightly ionized gas, net current 
of zero, and constant electron temperature. If the electron temperature is not constant, 
then E is proportional to the gradient of electron pressure and not to the electron density. 
Thus, Equation (2) is not valid in the shock region. 

In view of these difficulties, a new method for handling plasmas with DSMC is proposed. 
The concept of ambipolar diffusion is used explicitly in the modeling to determine charged 
particle motion and the electric field. Ambipolar diffusion results when the lighter electrons 
tend to diffuse faster than the ions for flows involving a mass density gradient. A charge 
separation and resulting electric field are produced. The electric field retards the electron 
diffusion while enhancing the diffusion of the positive ions. The modeling proceeds as follows. 

During the movement routine, the velocities of a charged particle in a given cell are 
determined from 

m s dvs/dt = q s E (3) 

The average velocity of each charged particle in a local region is written according to this 
equation. Then, the average velocity of ions and of electrons is determined by summing 
over the charged particles in that region. The resulting equations can be solved for the 
local electric field when the requirements of a net current of zero and charge neutrality 
are imposed. Details of this approach are given in the appendix. The present simulation 
accounts fully for both density and temperature gradient effects on the electric field. 

The electric field calculations are performed in the domain of a supercell, which consists 
of several adjacent computational cells (typically about 10). The supercell must contain 
several hundred particles in order to obtain a reasonable sample of charged particles. In 
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addition, the electric field result is time averaged before it is used in the calculation of 
charged particle motion. Because of the disparity between the collision frequencies of the 
electrons and the heavy particles, each is moved based on its own time step. This yields 
At e A th and electrons move more often than the heavy particles. The electric field, E, 

is calculated after each electron movement. The result is smoothed over the heavy particle 
time step, A th, and used for the next A th in the calculations of charged particle velocity. 
Once the flow has reached steady state, the electric field is calculated by averaging. 

Charge neutrality must be enforced in the simulation because the Debye length is much 
smaller than any other characteristic length. In the present approach, the electrons and 
ions are not tied together and charge neutrality is not automatically ensured. The following 
procedure is used if, after the movement routine is exercised, the charge neutrality condition 
is not satisfied. Randomly selected electrons in supercells with excess electrons are moved 
to randomly chosen locations in neighboring supercells which have a deficiency of electrons 
until the number of electrons and ions in each supercell is equal. 

Electron Impact Ionization 

The flow fields considered in this study are sufficiently energetic that a significant 
contribution to the total ionization is from electron impact ionization reactions. The 
inclusion of these reactions in the DSMC program presents certain problems. The reactions 
have traditionally been handled as single step reactions, but the available rate constants for 
these reactions are not of a form which can be used directly with the DSMC methodology. 

The continuum rate coefficients k(T) for chemical reactions are specified by 

k(T) = aT h exp(-E a /(kT )) (4) 

In the DSMC method, these rate coefficients are used to determine collisional energy 
dependent steric factors. The steric factor, Pr, is the ratio of the reaction cross-section 



to the total cross-section that results in the above rate coefficient. It is proportional to 


Pr <x (1 - E«/E c )< +t+1/2 , E e > E a 


(5) 


where ( is a measure of the vibrational and rotational internal degrees of freedom which 
may contribute to the reaction. Because Pr is zero for E c < E a , the representation is valid 
as long as 

( + b+ 1/2 >0 (6) 

The values of b given by Park and Menees[5] for the electron impact ionization impact 
reactions 

O + e >Q + +e + e (7) 

Af + e s- Af + + e + e (8) 

are —3.9 and —3.82 respectively. Because the reactants are monatomic gases, ( = 0 and the 
criterion indicated by Eqn.(6) is violated. The method that was employed to produce the 
rates traditionally used for these reactions involves determining an average flow temperature 
and replacing T b in the rate coefficient by its value at that temperature. This, in effect, 
changes k(T) to 

k(T) = Aexp(-E a /(kT )) , A = a(T avg ) h (9) 

In the DSMC simulations this is not accurate in the immediate vicinity of the shock 
because the temperature is changing rapidly and is considerably different from the average 
temperature used in the rate calculation. 

An alternate method is proposed which uses the assumption that these two reactions 
proceed via a two-step chain involving excitation followed by ionization from the excited 
state. Thus, the ionization of atomic nitrogen proceeds by 


Af A e > Af* + e 


( 10 ) 


Af* + e — > Af + + e + e 


( 11 ) 
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with similar reactions for atomic oxygen. In his paper on ionization in air behind high 
speed shock waves, Wilson[9] asserts that the rate limiting step in this ionization process 
is the excitation of nitrogen atoms to the 3 s 4 P state and oxygen atoms to the 3 s 5 S state. 
Experimental reaction cross-sections are available for these excitation processes from Stone 
and Zipf [1 0- 11]- These cross-sections are used directly in the program when determining 
the probability of excitation after a collision. 

The data show a much larger excitation cross-section for nitrogen than for oxygen for 
these two states. This suggests a higher rate of ionization for nitrogen, but the available data 
do not indicate that this is the case. An investigation of the processes involved reveals that 
the nitrogen 3 s 4 P state has a short radiative lifetime while the oxygen 3 s 5 S state is stable. 
Thus, the nitrogen may radiate from the excited state before it has a chance to ionize. To 
determine if the nitrogen atom radiates, a process similar to that used for rotational and 
vibrational transitions is employed. A relaxation collision number, R } is determined from 
the product of the radiative lifetime and the average collision frequency for nitrogen atoms 
at that point in the flow. The quantity 1/R gives the probability of transition for a single 
collision. For extremely rarefied flows the value for R may be less than 1. This indicates that 
the radiative lifetime is less than the average time between collisions and the nitrogen atom 
has very little chance of ionizing by this mechanism. Therefore, this reaction is bypassed 
if the relaxation collision number is less than 1. When the DSMC results using the one- 
step reaction rates for the same flow conditions are examined, it is found that the electron 
impact ionization of nitrogen occurs very rarely. Thus, no major discrepancy between the 
two methods is introduced by bypassing this reaction for very low density flowfields. 

For ionization from the highly excited state, the quantum defect method[12] and the 
experimentally determined cross-section equation of Lotz[13] are used. For a single electron 
in the subshell and an impact electron energy near threshold, Lotz gives 

a r = a(E/P 1 - 1)(1 - b)/Pl (12) 
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The values for the constants in the equation are determined in accordance with the 
assumption that the excited states are almost hydrogenic (a = 4. X 10 -14 cm/sec, b = .6). 
This is valid for all highly excited states of atoms. The quantum defect method requires 
the formula to be multiplied by a factor of n 4 where n is the effective principal quantum 
number. 

n = (Ry/(Eoo - Pi )) 1/2 (13) 


Thus, 


Or = n 1 a(E/P 1 — 1)(1 — b)/P? 


( 14 ) 


Results 

Results are presented for the 10 km/sec standing shock wave in air at .1 torr. 
This flowfield is of interest because it represents the conditions of an AVCO shock tube 
measurement [14], It is also representative of the flows which might be experienced by the 
AFE. In the plots, the flow is from negative to positive x with zero at the approximate 
shock center. For these conditions, the flowfield composition through the shock and the 
translational temperature profile are given in figures 1 and 2. 

One of the major purposes of the paper is to show that the DSMC represents a self- 
consistent framework and does not require borrowing any results from continuum theory in 
order to describe slightly ionized gases. Thus, when the present results are compared with 
those of Reference 6, the objective is not to pass judgement on the accuracy or suitability of 
the Langmuir- Tonks formala. We know it is inadequate when electron temperature gradients 
are important. Rather, it is to compare with the only other particle model that was used 
to calculate the flow under consideration. 

The present method for the plasma calculations has a significant effect on the electron 
temperature through the shock (figure 3). The peak temperature is lower and occurs further 
upstream of the shock. Also, energetic electrons are present through a much larger portion 
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of the flow. The sharp peak in electron temperature which is predicted with Bird’s method 
is the result of limiting the electron mobility by requiring the charged particles to move only 
in electron-ion pairs. The broader distribution of energetic electrons seems more physically 
reasonable. The shift in the peak demonstrates the tendency of the energetic electrons to 
diffuse to the lower density portions of the flow. 

The electric field for the two cases is shown in figure 4. The present method shows an 
electric field which is larger upstream of the shock than that which was predicted by the 
traditional method. This result agrees with the altered electron temperature distribution. 
The field is strongest where the electron velocity is largest. The difference in the magnitude 
of the electric field may result from the use of the continuum and constant electron 
temperature assumptions of the Langmuir- Tonks equation used in Bird’s method. 

The electron concentration is somewhat less for the present method than for the Bird’s 
method (figure 5). This is probably because the peak electron temperature is lower and 
occurs further upstream where the heavy particle temperature is lower. As a result there is 
less energy available for reactions involving electron impact in this portion of the flow. 

Implementing the two-step electron impact reaction model has some effect on the results 
for electron temperature in the immediate vicinity of the shock (figure 6). The two-step 
modeling allows more accurate matching of the experimental cross-sections and reaction 
rates in the high temperature gradient region. The electron density (figure 7) is a bit 
smaller at the center of the shock but increases in the region behind the shock. As can be 
seen from the concentration of nitrogen and oxygen ions (figures 8 and 9) the ionization of 
the atoms is increased in the region approximately .25 to .5 cm behind the shock when the 
two-step modeling technique is used. 


Conclusions 

Some changes to the DSMC method for modeling charged particles and electric 
field effects have been introduced. These changes allow a more accurate representation 
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of ambipolar diffusion effects and electron impact ionization reactions. The results are 
noticeably different for electron density and electron temperature in the flow. It is 
recommended that this method be used in future simulations when the determination of 
these properties is essential. 
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Appendix: Calculation of Electric Field 


The instantaneous electric field in a given electric field cell (or supercell) is determined 
from the consideration that the net ion current is equal to the net electron current. The 
resulting quantity, like many others of interest such as heat flux, stress, etc., is given in 
terms of the instantaneous particle velocities. The mean electric field at a given point in 
the flowfield is determined by averaging over a large sample. 

Replacing the electric field by its average value over the time step At, the average 
velocity of a charged particle during At follows from equation (3) as 

v s = v Sj0 + (q s EAt/2m s ) (15) 


where v Sj o is the particle velocity at the beginning of the time step, m s is the particle mass 
and q s is the particle charge. The average ion velocity for that time step is 








(16) 


where a refers to an ion species and N a is the number of particles of that species in the 
supercell. The average electron velocity during the time step is 


V, = - ( eE/m,){At/2 ) 


(U) 


where e is the charge of the electron. Because charge neutrality is required in the supercell, 
N e = 22 N a , assuming singly ionized ions. The electric field is calculated by setting V = V e 
and solving the resulting equation. The result can be written as 


eEAt/2 


E^e,p - 22^,0 ]/N e 
22 1 /m a + l/m e 







( 18 ) 


The above formula for the instantaneous electric field is averaged over many samples 
to give the mean electric field. 
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Captions 


Fig. 1 Composition through 10 km/s shock in air at .1 torr 

Fig. 2 Translational temperature through the shock 

Fig. 3 Electron temperature through the shock 

Fig. 4 Electric field through the shock 

Fig. 5 Electron concentration through the shock 

Fig. 6 Electron temperature through the shock 

Fig. 7 Electron concentration through the shock 

Fig. 8 N+ concentration through the shock 

Fig. 9 0 + concentration through the shock 
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